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This article describes the Schwarzschild orbit superposition method. It is the state-of-the-art dynamical modelling tool for 
early-type galaxies. Tests with analytic models show that masses and orbital anisotropies of not too face-on galaxies can 
be recovered with about 15 percent accuracy from typical observational data. Applying Schwarzschild models to a sample 
of Coma galaxies their dark matter halos were found to be 13 times denser than those of spirals with the same stellar mass. 
Since denser halos assembled earlier, this result indicates that the formation redshift 1 + Zf orm of ellipticals is about two 
times higher than of spirals. Roughly half of the sample galaxies have halo assembly redshifts in agreement with their 
stellar-population ages. Galaxies where stars appear younger than the halos show strong phase-space density gradients in 
their orbital structure, indicative for dissipational evolution and possibly connected with secondary star-formation after the 
main halo assembly epoch. The importance of considering dark-matter in dynamical models aimed to measure black-hole 
masses is briefly discussed. 
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1 Introduction 

Early-type galaxies are characterised by an overall smooth 
and featureless spheroidal morphology and a dynamically 
hot system of stellar orbits. This is thought to be the result 
of a dynamically violent assembly process. However, when 
and how exactly these galaxies have formed is still not well 
known. Their mostly old and a-enhanced stellar populations 
imply a relatively short star-formation period in the distant 
past. The epoch when those stars assembled to form the 
early-types seen today cannot be directly deduced from stel- 
lar population ages. For example, the stars might be born at 
high redshift, in progenitors that only recently merged into 
spheroidal galaxies. If these mergers are mostly collision- 
less (without significant amounts of gas and star- formation), 
then stellar-population ages do not change, and the main as- 
sembly epoch is delayed with respect to the star-formation 
epoch. In contrast, a spheroidal galaxy formed early in the 
universe could have grown a disky stellar subcomponent re- 
cently, e.g. triggered by a gas-rich minor merger. The main 
halo assembly would then precede star formation (along the 
disk). 

Pure dark-matter, collisionless iV-body simulations pre- 
dict a close relationship between the main assembly redshift 
of dark matter halos (when, say, half of the mass had been 
assembled) and their average density. The infall of baryons 
might change the dark matter density as it enforces an ex- 
tra gravitational pull on the halo. Still, measuring the dark 
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matter density in elliptical galaxy halos is a valuable tool 
to gain information about their assembly epoch. In compar- 
ison with stellar-population ages it also gives indirect in- 
formation about the evolution of these systems. According 
to the above considerations one can expect that galaxies in 
which stars appear younger than the halos have likely ex- 
perienced some secondary star-formation episodes. These 
systems have evolved dissipatively. In contrast, when stars 
appear older than the halo, this might indicate a more gas- 
less evolution (e.g. by gas-free or collisionless or dry, re- 
spectively, mergers). 

Since the system of stars in galaxies is collisionless it 
preserves some information about how the stars have as- 
sembled. A dynamically violent formation - characterised 
by violent relaxation in phase-space - is supposed to result 
in a highly mixed orbit distribution with strong phase-space 
density gradients likely being washed out. In contrast, dis- 
sipational evolution (e.g. through gas-rich mergers) likely 
results in disky subsystems with high phase-space density 
peaks on near-circular orbits. Consequently, the analysis of 
the orbital structure provides additional information about 
the assembly mechanism of early-type stellar systems. 

In the last years we have collected photometric and kine- 
matic observations for a sample of 19 early-type galaxies in 
the Coma cluster: 2 central cD galaxies, ten giant ellipticals 
and seven SO or E/S0 galaxies, respectively, with —18.8 > 
M B > -22.26 (Mehlert et al. 2000; Wegner et al. 2002; 
Corsini et al. 2008; Thomas et al. 2009a). For all galaxies 
ground-based and HST photometry are available. Long-slit 
stellar absorption-line spectra have been taken along at least 
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the major and minor axes. In many cases additional data 
along other position angles has been collected as well. The 
spectra extend to 1 — 4 r e ff, entering the regions where dark 
matter becomes noticeable. 

By means of dynamical models we have measured the 
dark matter content and orbital structure of these galaxies. 
The results are summarised below and implications for the 
formation process of early-type galaxies are discussed. 

2 The Schwarzschild method 

A complete description of a stellar system is provided by 
its phase-space distribution function /, i.e. the density of 
stars in 6-dim phase-space. Unlike for a collisional gas the 
distribution function of collisionless dynamical systems like 
galaxies is not known in advance. However, for steady-state 
objects Jeans theorem ensures that the phase-space density 
is constant along individual orbits. Orbits, in turn, are iden- 
tified by the integrals of motion they respect. Thus, for sys- 
tems in a steady state the distribution function (DF) reads 

f = f(h,...,I n ), (1) 

where I\,. . . ,I n are integrals of motion (e.g. Binney & 
Tremaine 1987). The simplest symmetry assumption con- 
sistent with the flattening and rotation of elliptical galaxies 
is that they are axisymmetric. Orbits in typical axisymmetric 
galaxy potentials respect three integrals of motion: energy 
E, angular momentum along the symmetry axis L z (the ro- 
tation axis being parallel to the z-axis) and the so-called 
third integral I3 (e.g. Contopoulos 1963). The last integral 
is usually not known explicitly and the distribution function 
can thus not be written in terms of elementary functions. 

Schwarzschild (1979) introduced an orbit-superposition 
technique (now called Schwarzschild method) to construct 
collisionless DFs. In brief, one assumes a trial potential for 
a given galaxy and composes a library of several thousand 
orbits. The light distribution and projected kinematics of 
each orbit are stored. Then, a galaxy model is constructed 
as the superposition of all orbits. The unknown weight or 
total amount of light, respectively, of each orbit is chosen 
to match the model as good as possible with the given con- 
straints (see below). This method corresponds to the approx- 
imation / w ^ li fi, where the fa are single-orbit DFs (Van- 
dervoort 1984; Thomas et al. 2004). The accuracy of the 
method only depends on the density of the orbit grid. 

The Schwarzschild method provides very general dy- 
namical models. In contrast to Jeans models there is no ne- 
cessity for any a priori restriction upon the orbital structure. 
Moreover, Schwarzschild models are easily constructed to 
be everywhere positive in phase space (i.e. physically mean- 
ingful) - a property that is not guaranteed in Jeans models. 
The main challenge is to ensure that the orbit library is rep- 
resentative for all orbital shapes. 

In axisymmetric potentials it is straight forward to sam- 
ple the energy E and the angular momentum L z . The issue 
related to our ignorance about I3 is overcome by a system- 
atic sampling of orbital initial conditions, which implicitly 



guarantees the inclusion of all orbital shapes. More specif- 
ically, orbits with given E and L z but different I3 follow 
distinct invariant curves in appropriate surfaces of section 
(SOS). Launching orbits at constant E and L z from vari- 
ous initial starting points on a grid in such SOSs ensures a 
representative sampling of the unknown ^3 (Thomas et al. 
2004). 

In practice, one goes through the following steps to con- 
struct a Schwarzschild model of a galaxy: 

- The photometry of the galaxy is deprojected to obtain 
the 3d light distribution v. 

- A trial mass distribution p is set up, combining the vari- 
ous mass components, e.g. 

p = T x v + pdm + Mbh x S(r). (2) 

Here, T is the stellar mass-to-light ratio (typically as- 
sumed to be radially constant in ellipticals), pdm is the 
dark matter density (see below) and Mbh represents a 
supermassive central black hole. Given p, the gravita- 
tional potential follows by solving Poisson's equation 
and an orbit library can be assembled. 

- The orbits are superposed and the orbital weights that 
best match with the observations are determined. In our 
implementation this is done by maximising 



S — a\ 



(3) 



where \ 2 quantifies deviations between observed and 
modelled kinematics (the observed light profile is used 
as a boundary condition in the maximisation). The func- 
tion 

S=- f f\nfd 3 rd 3 v (4) 

is the Boltzmann entropy of the orbit distribution and 
a is a regularisation parameter (Richstone & Tremaine 
1988; see also below). 
- The parameters determining the potential (e.g. T, halo 
parameters, Mbh) are systematically varied and the fi- 
nal best-fit is obtained from a \ 2 analysis. 
Various implementations of this method exist for spheri- 
cal (Romanowsky et al. 2003), axisymmetric (Cretton et al. 
1999; Cappellari et al. 2006; Chaname, Kleyna & Van der 
Marel 2008; Valluri, Merritt & Emsellem 2004; Gebhardt 
et al. 2003; Thomas et al. 2004) as-well as triaxial poten- 
tials (van den Bosch et al. 2008). Gebhardt et al. (2003) 
measured supermassive black-holes in the centres of galax- 
ies (ignoring any contribution of dark matter). Cappellari et 
al. (2006) modelled a subsample of the SAURON galaxies, 
again ignoring dark matter because their data cover only the 
inner regions < r c s- 

Beyond r e s dark matter becomes important in early- 
type galaxies. The Coma sample is by now the only larger 
sample of generic (e.g. flattened and rotating) ellipticals that 
has been modelled with Schwarzschild's method including 
dark matter. Previous attempts to measure the dark mat- 
ter content of early-types via stellar dynamics focussed on 
round and non-rotating systems and assumed spherical sym- 
metry (Kronawitter et al. 2000; Gerhard et al. 2001 ; Magor- 
rian & Ballantyne 2001). 
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Fig. 1 Dark matter density (pdm) (averaged within 2 r en 
versus stellar mass M*. Red dots: Coma early-type galax- 
ies; open circles: round and non-rotating galaxies from Ger- 
hard et al (2001); small dots: semi-analytic galaxy forma- 
tion models (De Lucia & Blaizot 2007). Blue lines: spiral 
galaxies. 



In the Coma galaxies we probed for two different para- 
metric halo profiles: (1) logarithmic halos with a constant- 
density core and asymptotically constant circular velocity 
and (2) so-called NFW-profiles which are good fits to dark 
matter halos in cosmological A^-body simulations (Navarro, 
Frenk & White 1996). The majority of Coma galaxies are 
better fit with logarithmic halos, but the significance over 
NFW halo profiles is marginal (Thomas et al. 2007b). 

3 Tests of the method 

A unique feature of our Schwarzschild models is the incor- 
poration of orbital phase-space volumes V (Thomas et al. 
2004). Together with the total amount of light w on the or- 
bit it allows to calculate the phase-space density / = w/V 
on each orbit. This offers several new applications. Firstly, 
one can project analytic phase-space DFs /analytic through 
orbit libraries. In this case, the orbital weights are not deter- 
mined from a fit to photometric and kinematic constraints, 
but directly from the analytic DF, i.e. the weight Wi of or- 
bit i reads Wi = /analytic x Vi. The projected kinematics, 
spatial density profile, and intrinsic velocity dispersion pro- 
files of the so-constructed orbit superposition can be com- 
pared with direct phase-space integrations of /analytic- Both 
methods agree very well (Thomas et al. 2004), confirming 
that the orbit sampling is representative. 



In addition, mock data from analytic model galaxies can 
be used to measure how accurate the orbital DF can be re- 
constructed. In this Monte-Carlo approach photometric and 
kinematic data with spatial resolution and coverage simi- 
lar to real data are simulated and modelled exactly as real 
galaxies. For the Coma data we have shown that with the ap- 
propriate choice of the regularisation parameter a the orbital 
structure and mass distribution of not too face-on galax- 
ies can be recovered with an accuracy of about 15 percent 
(Thomas et al. 2005). Similar tests have been presented in 
Cretton et al. (1999), Krajnovic et al. (2005) and Siopis et 
al. (2009). 

The input models for the above tests had the same sym- 
metry as the dynamical models (axial symmetry). Isophotal 
twists and kinematic misalignment suggest that at least the 
most massive early-type galaxies are slightly triaxial. In or- 
der to explore the systematic errors arising from too restric- 
tive symmetry assumptions, we applied our models also to 
mock data sets created from collisionless A^-body binary 
disk merger remnants. These remnants are strongly triaxial 
in the central, box-orbit dominated regions (Jesseit, Naab 
& Burkert 2005). While the central mass measured with ax- 
isymmetric models then underestimates the true mass on av- 
erage (depending on projection angle and box-orbit content) 
the enclosed mass within r c g is still recovered mostly with 
better than 20 percent accuracy unless for highly flattened, 
face-on systems (Thomas et al. 2007a). 

4 The dark matter density and assembly 
epoch of early-type galaxies 

Fig. Q] shows dark matter densities (pdm) as a function of 
stellar mass. Stellar masses A/* = T x L of early-types are 
derived from the best-fit (dynamical) stellar mass-to-light 
ratio T and the observed total luminosity L. Dark matter 
densities (pdm) are averaged within 2 r c ff . 

The dark matter densities of Coma ellipticals (from ax- 
isymmetric modelling) and round and non-rotating galax- 
ies (from spherical models) match well. With increasing 
stellar mass dark matter densities tend to decrease. Spiral 
galaxy dark matter densities are also shown in Fig. [T] (Per- 
sic, Salucci & Stel 1996a,b; Kormendy & Freeman 2004, 
respectively). Stellar masses for spirals have been derived 
using the Tully-Fisher and stellar-mass Tully-Fisher rela- 
tions of Bell & De Jong (2001) (cf. Thomas et al. 2009a 
for details). As in early-types, dark matter densities in spi- 
rals decrease with increasing stellar mass. Most important, 
the dark matter in ellipticals is about 13 times denser than 
in spirals of the same stellar mass (compared at the same 
_B-band luminosity, early-type galaxies have 7 times denser 
halos than late-types). 

Finally, Fig.Q]also includes semi-analytic galaxy forma- 
tion models of De Lucia & Blaizot (2007). They are in good 
agreement with the observations in Coma galaxies. This is 
surprising because the A^-body cosmological simulation un- 
derlying the semi-analytic models was performed without 
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Fig. 2 Halo assembly redshifts (y-axis) versus star- 
formation redshift (x-axis) for Coma early-type galaxies. 
The one-to-one relation is indicated by the dotted line. 
The phase-space DFs of the two galaxies NGC4827 and 
NGC493 1 are shown in Fig. @] 



baryon dynamics. Thus, either the net effect of baryons on 
dark matter around early-types in the given mass interval 
vanishes or there is some discrepancy between the semi- 
analytic models and the observations. Note that baryons can 
also cause a halo expansion (e.g. by dynamical friction). 

The light distribution in late-type galaxies is less cen- 
trally concentrated than in early-types. Then, even if the net 
effect of baryons on elliptical galaxy halos would be negli- 
gible, the stronger gravitational pull in early-type galaxy ha- 
los could still contribute to the observed overdensity of dark 
matter in ellipticals relative to spirals. Within the adiabatic 
contraction approximation this can at most explain a factor 
of two between average halo densities around ellipticals and 
spirals, respectively (Thomas et al. 2009a). The remaining 
over-density indicates that ellipticals have assembled earlier 
- at a time where the universe was denser. 

The simplest analytic models, as well as pure dark mat- 
ter cosmological iV-body simulations predict a scaling of 
the average dark matter density (pdm) with halo-assembly 
redshift z form of (/9dm) oc (1 + Zform) 3 (e.g. Gerhard et al. 
2001). Accordingly, the density contrast between elliptical 
and spiral galaxy halos translates into 1 + Zf orm being about 
two times higher for ellipticals than for spirals. Absolute 
assembly redshifts can only be estimated with an additional 
assumption upon the formation redshift of spirals. Supposed 
that spirals assemble typically at Zform ~ 1 ( at higher red- 
shifts spirals become rare, e.g. Conselice et al. 2005) then 



Coma early-types have assembled around Zf orm ra 2 — 3. 
These dark-matter based formation redshifts are shown in 
Fig. |2] against stellar-population ages. In about one third of 
Coma galaxies the stars appear to be younger than the halos 
but many galaxies are consistent with equal assembly and 
star- formation redshifts. 

5 The orbital structure of early-type galaxies 

As outlined above a galaxy in which the stars appear to 
be younger than the halo is a candidate dissipative system, 
while objects in which stars are equally old or older than the 
halo have likely formed monolithically or through mostly 
collisionless mergers. In dissipative systems one would ex- 
pect a high phase-space density on near-circular disk orbits 
while the phase-space DF of systems which have undergone 
violent relaxation should lack strong phase-space density 
gradients. 

Classically, the orbital structure is measured in terms of 
anisotropy parameters, i.e. ratios of intrinsic velocity dis- 
persions. Let / denote the phase-space DF of a galaxy, then 
the intrinsic dispersions <7ij read 

f(vi -vt){vj -vj)d 3 v, (5) 



9 



Vj 



/«, 



a\ 



(6) 



In the following we assume i,j £ {R, </>, z}, where z is 
the symmetry axis (cylindrical coordinates). The unordered 
kinetic energy along coordinate direction i is 



n,, ; = 



P4 dV 



(7) 



and the anisotropy in the velocity dispersions can be quan- 
tified by j3 = 1 - H zz /n RR and 7 = 1 - U^/Hrr (Cap- 
pellari et al. 2007). An isotropic system has j3 — 7 = 0. 

Fig.[3]shows (3 and 7 versus intrinsic ellipticity e. In ob- 
served galaxies > and, on average, 7 sa 0, but with 
significant scatter. The range of anisotropies found in Coma 
early-types is similar as in SAURON galaxies but the Coma 
galaxies do not exhibit a trend of increasing (3 with increas- 
ing e as observed by Cappellari et al (2007). Neither the 
inclusion of dark matter in the Coma models nor regular- 
isation can explain this difference (Thomas et al. 2009b). 
Instead it most likely reflects differences in the selections of 
the two samples (Thomas et al. 2009b). 

In order to clarify the relationship between the classi- 
cal anisotropy and the intrinsic orbital structure two proto- 
typical orbital compositions are shown by the lines in Fig. [3] 
(1) The case of flattening by rotation, i.e. by extra-light on 
near-circular orbits. Circular orbits mostly contribute to the 
kinetic energy in <\> direction and 7 becomes negative. The 
more, the flatter the galaxy is (dashed lines). (2) The dot- 
ted lines show toy models with maximum entropjQ With- 
out any other conditions the maximisation of Boltzmann's 



1 For both toy models 75 percent of the light was distributed on prograde 
orbits; cf. Thomas et al. (2009b) for details. 
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Fig. 3 Anisotropy j3 and 7 versus intrinsic ellipticity e 
for observed galaxies (left-hand panels) and A^-body binary 
disk merger remnants (right-hand panels). Dotted lines: 
maximum-entropy toy-models, dashed lines: 21 models (de- 
tails in the text). The phase-space DFs of the two galaxies 
NGC4827 and NGC4931 are shown in Fig.S 



entropy yields a fiat DF / = const. This is altered as soon 
as the luminosity distribution is used as a boundary con- 
dition. Still, maximising the entropy yields, in a sense, the 
smoothest phase-space DF compatible with a given density 
profile. It turns out that the classical notion of flattening by 
anisotropy (i.e. a suppression of vertical versus horizontal 
energy, with isotropy between R and <f>) is closely related 
to the maximisation of entropy. Accordingly, systems along 
the maximum-entropy line, while having varying anisotropy 
j3, are similar in that their phase-space distribution functions 
are smooth. 

The phase-space DFs of two Coma galaxies shown in 
Fig. 2] (both galaxies are flagged by the arrows in Figs. [2] 
and 0. The two galaxies have similar j3 and 7 but differ- 
ent intrinsic flattening: NGC4827 is an E2 early-type while 
NGC493 1 is highly flattened. In addition, in NGC4827 stars 
are about as old as the halo while the stars in NGC493 1 ap- 
pear younger than its halo. 

In the upper panels of Fig. [4] spheroidal orbits are high- 
lighted. Disk orbits are highlighted in the lower panels. We 
identify spheroidal orbits by |^| max > 70°, where d is the 
angle between a point on the orbit and the equatorial plane 
of the galaxy. So-defined spheroidal orbits come close to the 
galaxy's minor-axis and contrast disk orbits, which are in- 
stead selected by |z| ma x < r e ff/4, where z is the vertical 
height of a point on the orbit with respect to the equatorial 
plane. 

Fig.@]shows that phase-space densities on spheroidal or- 
bits are similar in both galaxies. Moreover, in NGC4827 the 
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Fig. 4 Phase-space distribution function / for the two 
galaxies NGC4827 (left-hand panels) and NGC4931 (right- 
hand panels). Each dot represents the phase-space density 
on a single orbit (in solar masses per cubed kpc and cubed 
km/s) plotted against the mean orbital radius. In the top pan- 
els spheroidal orbits are highlighted; in the bottom panel 
disk orbits are highlighted (details in the text). Only pro- 
grade orbits are shown. 



phase-space densities of disk orbits are comparable to those 
of spheroidal orbits. Yet, the stellar density on disk orbits in 
NGC493 1 is up to two orders of magnitude higher than on 
spheroidal orbits. These findings support the interpretation 
that galaxies left to the one-to-one relation in Fig. [2] (stars 
younger than halo) are dissipative systems, while galaxies 
close to the one-to-one relation or right of it are systems 
which formed in a dynamically violent process and have 
smoother phase-space DFs. A more detailed investigation 
of phase-space densities will be presented in Thomas et al. 
(in preparation). 



The right-hand panels of Fig. |3]display dynamical mod- 
els of collisionless A^-body binary disk merger remnants. 
Mock data sets with similar spatial resolution and cover- 
age as for Coma galaxies were constructed from projections 
along the three principal axes of six remnants selected rep- 
resentatively from the A^-body sample of Naab & Burkert 
(2003). While the distribution of /3 versus e is like in ob- 
served galaxies, none of the analysed merger remnants has 
7 < 0. In view of the above discussion this is plausible since 
galaxies with 7 < are likely dissipative systems (large en- 
ergy in <p caused by extra-light on circular orbits), while the 
analysed A^-body merger simulations are gas-free. 
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6 Summary 

The Schwarzschild technique is the state-of-the-art tool to 
model early-type galaxies. For typical observational data 
masses and orbital anisotropies can be recovered with about 
15 percent accuracy (for not too face-on systems). We have 
applied Schwarzschild models to 19 Coma early-type galax- 
ies to measure the dark matter content and distribution of 
stellar orbits. By today, it is the largest sample of generic 
early-type galaxies with dynamical models including dark 
matter. Dark matter densities in early-types are larger than 
in spirals of the same luminosity or mass. Extra gravita- 
tional pulling by the more concentrated baryons in ellipti- 
cals is not sufficient to explain this over-density. Instead, the 
formation redshift l+Zf or m of early-types is about two times 
higher than of spirals. Under the assumption Zf orm w 1 for 
spirals the Coma early-types have formed around Zf orm w 
2 — 3. Observed dark matter densities are in good agreement 
with recent semi-analytic galaxy formation models. In about 
half of the sample galaxies halo assembly redshifts match 
with stellar population ages. In galaxies where stars appear 
younger than the halo, the orbit distribution indicates dis- 
sipational evolution (i.e. strong phase-space density peaks 
on near-circular disk orbits). This suggest that these galax- 
ies had some secondary star-formation after the main halo 
assembly epoch. 

When modelling the very central regions of early-types 
to study their supermassive central black-holes it is impor- 
tant to include dark matter in the models, even if the central 
parts itself are not dominated by dark matter. Neglect of the 
halo can result in an overestimation of the stellar mass-to- 
light ratio which subsequently leads to an overestimation of 
the central stellar mass. The latter is degenerate with the 
black-hole mass and neglect of dark matter might then yield 
a too small black-hole. This has been illustrated for M87, 
where the black-hole mass more than doubled after includ- 
ing dark matter in the models (Gebhardt & Thomas 2009). 
The effect is supposed to be strongest for the most mas- 
sive galaxies since their light profiles are shallow. Further 
galaxies have to be modelled in order to establish if the in- 
clusion of dark matter in black-hole models can reduce the 
discrepancy between the s=s 10 9 M Q solar mass black-holes 
in the most massive nearby galaxies and the w 10 10 M & 
solar mass black-holes in high-redshift quasars. 
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